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Abstract. The classical two-dimensional discrete frustrated (j>^ model is studied by Monte Carlo simula- 
tions. The correlation function is obtained for two values of a parameter d that determines the frustration 
in the model. The ground state is a ferro-phase for d = —0.35 and a commensurate phase with period 
A*' = 6 for d = —0.45. Mean field predicts that at higher temperature the system enters a para-phase 
via an incommensurate state, in both cases. Monte Carlo data for d — —0.45 show two phase transitions 
with a floating-incommensurate phase between them. The phase transition at higher temperature is of the 
Kosterlitz-Thouless type. Analysis of the data for d = —0.35 shows only a single phase transition between 
the floating-fluid phase and the ferro-phase within the numerical error. 



PACS. 64.70.Rh Commensurate-incommensurate transitions 
vibrations and displacive phase transitions 



63.70.-|-h Statistical mechanics of lattice 



1 Introduction 

One of the microscopic origins of the spatially modulated 
structures is given by the competing interactions between 
the particles. The Frenkel-Kontorova model, the axial next- 
nearest-neighbour Ising (ANNNI) model, and the discrete 
frustrated (DIFFOUR) model provide basic microscopic 
Hamiltonians for the systems of this class. These models 
have been studied extensively, specially the first two PH]. 
Competing interactions lead to the appearance of incom- 
mensurate and commensurate regions in the phase dia- 
grams of these models. The order parameter in these sys- 
tems is a vector that determines the amplitude and the 
phase of modulation. For the incommensurate state, the 
free energy is degenerate in the phase of modulation. The 
two-dimensional case is of special interest here, since there 
is no long-range order in such planar systems if the order 
parameter has more than one component |3| ■ A nonmono- 
tonic algebraic decay of correlations appears instead of the 
usual incommensurate state 

The correlation function in two-dimensional systems 
with competing interactions has at least three possible 
types of behaviour, as found in the spin models 2 . The 
disordered phase shows an exponential monotonic or mod- 
ulated {q 0) fall-off of the correlation function: 

C{r) =< S{0)S{r) cxp(-r/0 cos{qr). (1) 

The case g 7^ is usually called floating-fluid phase. An 
algebraical decay of correlations at r — > 00 is a property of 
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the floating-incommensurate (FIC) or floating-solid phase: 



C(r) ~ r ^' cos{qr + ^p) . 



(2) 



Finally, a commensurate phase with locked-in wavevector 
(Jo is not destroyed by the fluctuations in two-dimensional 
case. Therefore the correlation function in this case takes 
the form 

C{r)^cos{qor). (3) 

The commensurate case includes the ferro-phase also (go = 
0). 

One should note that in the two-dimensional case the 
phase transition from the para-phase to the FIC-phase 
is believed to be related with the formation and disso- 
ciation of the dislocations or vortices This scenario 
of a phase transition can be described by the Kosterlitz- 
Thouless theory which is basic for the two-dimensional 
XY universality class The power law dependence 

of C(r) for this universality class appears everywhere be- 
low the critical temperature, down to T = or to another 
phase transition (to a commensurate state, for example). 
The critical index 77 is equal to 0.25 in the point of the 
phase transition (Tc) and decreases for decreasing temper- 
ature as ri{T) — 0.25 — C^/Tc — T, where C is a constant 
{rj ^ T for T ^ 0). The thermal dependence of the cor- 
relation length above the transition point is of the form 
^ ^ exp(6/ ^JT — Tc), where 5 is a constant. 

One of the well studied spin models with competing 
interactions is the ANNNI model. The ANNNI model is a 
spin lattice model with ±1 classical spins, ferromagnetic 
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nearest-neighbours interaction Ji in all directions and an- 
tiferromagnetic ncxt-nearest-neighbour interaction J2 in 
a single direction. The ANNNI model in two dimensions 
has been studied by numerous methods PE] ■ The ground 
state of this model at zero temperature is a ferromag- 
netic phase for — J2/J1 < 0.5 and a < 2 > phase (two 
spins up and two spins down ttii) for — J2/J1 > 0.5. Ob- 
viously, an FIC-phase can appear only at finite temper- 
ature between the paramagnetic phase and the < 2 > 
phase for —J2/J1 > 0.5. The behaviour of the model 
was studied by transfer-matrix calculations and finite-size- 
scaling analysis [Jj, Monte Carlo simulations [HlEl, dy- 
namical Monte Carlo method ^Jj, fermion approxima- 
tion ^] and its modification jj^lj the interface free en- 
ergy method 13 . Some of these results seem to indi- 
cate a Kosterlitz-Thouless 0] type of the phase transition 
from the paramagnetic to the FIC-phase and a Pokrovsky- 
Talapov ^3] type of phase transition from this phase to 
< 2 >. The latter scenario implies a square-root power law 
change of the wavevector with respect to temperature near 
the transition to the < 2 > phase. One should point out 
that the one-dimensional quantum ANNNI model studied 
in jTH] also reveals all phases and types of the phase transi- 
tions mentioned above. The approximations give a spread 
of the phase transition temperature over a rather wide 
region. On the other hand, other calculations predict an 
FIC-phase with a width comparable to that spread. The 
recent results on a cluster heat bath simulation § of sys- 
tems of about 64 x 128 atoms demonstrate that the FIC 
region is at least 2 times smaller than the previous pa- 
pers predicted 7 . Further, the nonequilibrium relaxation 
method JjQ^ for 6400 x 6400 systems does not predict any 
area for the FIC-phase at all, so that even the presence of 
this phase was questioned. 

In this paper we consider the classical two-dimensional 
DIFFOUR, which has continuous variables and compet- 
ing interactions in one direction. The DIFFOUR model 
was introduced in as a displacement translationally 
invariant model for the incommensurate phases. Instead 
of ±1 spins, this model considers double- well anharmonic 
2-4 oscillators at the nodes of a lattice. Similarly, there is 
a harmonic nearest-neighbour and next-nearest-neighbour 
coupling. The potential energy can be written as 

+ ^E("'-"^)'4- (4) 

Here cr.y = 1 for the nearest neighbours and cTy = 
otherwise, while a'^ — 1 for the next-nearest-ncighbours 
in the direction z (see Fig. la). The phase diagram for the 
analogous three-dimensional model is presented in |T5l . 
The system has commensurate and incommensurate re- 
gions in the phase diagram for d < i.e. with competing 
interactions. 

The ANNNI model can be considered as a limiting 
case of the DIFFOUR model (a +00). It corresponds 
to the order-disorder limit of the latter model. In the other 



limit (a +0) of the DIFFOUR model the ground state 
may be an incommensurate phase which is reached via a 
displacive transformation |19| . Furthermore, the greater 
freedom in the DIFFOUR model could also possibly lead 
to the existence of a broader the FIC-phase than in the 
ANNNI model. In view of the unclear situation in the 
latter we want to study the existence of FIC-phase in the 
DIFFOUR model. We shall do this using Monte Carlo 
simulations. The mean- field phase diagram can easily be 
determined and is presented in Fig. lb. 

Correlation functions of the model are analyzed for cer- 
tain parameters; their behaviour is always one of the three 
scenarios mentioned above (|1I2I3I) . An FIC-phase has been 
found. Indications of the Kosterlitz-Thouless scenario of 
the phase transition from the floating-fluid to the FIC- 
phase are discussed. The transition to the ferro-phase oc- 
curs via the floating-fluid phase. Surprisingly, the period 
of this phase increases as temperature decreases. 



2 Algorithm 

We study the two-dimensional DIFFOUR model using the 
Monte Carlo method. The difficulty of the simulation of 
the models with competing interactions is in the complex- 
ity of their configuration space which usually has many 
deep almost degenerate local minima. Then, it is difficult 
to overcome a potential barrier and to change the period 
of the modulation. This situation can lead to an unac- 
ceptably long computation time using conventional Monte 
Carlo schemes, particularly for large systems. Special al- 
gorithms should be used to achieve a reasonable result [S] 

IKl. 

In this paper, we use the "local-heating" algorithm 
which worked well for the three-dimensional DIFFOUR 
model JHl- The idea is to modify the potential energy 
in certain regions of the system to allow an easier switch 
between local minima. Actually, the temperature of the 
system is not changed but the constants of the potential 
energy are modified. However, for the largest part of the 
system its potential energy is still determined by Q . The 
regions with a modified potential energy are not included 
in the calculation of the physical quantities (order param- 
eter, correlation functions), while the usual Metropohs al- 
gorithm for Monte Carlo sampling is applied to the whole 
system. 

We use periodic boundary conditions and the following 
modification of the potential energy: 

i i i,j 

+ ^ E ^(«+j)/2("i - "jO^CTy - 

0, = {l + h{exp[-{x^^+zf)/q] 

+ exph((7V, -X,)' + {N, ~ z,f)/ll]))-\h > 0. (5) 

Nx and are the sizes of the lattice along the x and z- 
directions, respectively. The new coefficients 0i ((HJ mean 
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Fig. 1. a) Scheme of interactions in the two-dimensional DIFFOUR model, b) Mean-field phase diagram of the two-dimensional 
DIFFOUR model for the parameter a — 5 {t is the temperature hereafter). P - para-phase, F - ferro-phase, IC - incommensurate 
phase, C - commensurate phase. There are also narrow regions of the commensurate phases in the area of the incommensurate 
phase not shown in the picture. 



that their values in the vicinity of the lines Xi = and 
Zi = are smoothly decreased. Parameters h, I2 determine 
the height and the width of this decrease, respectively. 
The change of the parameters mimics a local heating of 
the corresponding regions. The idea of this procedure is 
to put these parts of the system in a para-phase and to 
destroy the local modulation phase. If the profile 6{x, z) is 
smooth enough, the phases of modulation at the two sides 
of the heated line are in an ideal case independent and 
can randomly fluctuate. The whole system may change 
the period of the modulation now. Use of " local heating" 
leads to faster thermalization of the system. 

Additionally we realised that it is useful to simulate 
the systems with a larger size along the x direction. The 
possible explanation is that once the period is changed at 
some X, it is relatively easy for this change to spread to 
the whole system. The probability for this initial change 
is proportional to the a;-size. 

We calculate the correlation function as follows: 



C(r) =<u(r')u(r' + r) > |, 



(6) 



Here <> |r' means an ensemble average. We consider two 
functions Cx = C{r = r^) (along x axis) and Cz = C(r = 
Tz) (along z axis) hereafter. The presence of the heated 
regions is taken into account as follows. We add data to 
the point r of correlation function lO if there is no " local 
heated" region between r and r' and add data to the point 
r — r' otherwise. If one of the points r or r' appears to be 
in the "local heated" region the data are not included in 
the calculation of the correlation function. This procedure 
reduces the influence of the finite-size effects. The data 
for the correlation function appear to be reliable even for 
distances larger than half the system size. 



We check whether the approach of the "local heat- 
ing" algorithm by changing the size of the system and by 
starting with various initial distributions of atom's dis- 
placements is reasonable. One should use large enough 
sizes of the system since it is difficult to distinguish expo- 
nential and power law decay of correlations. The situation 
becomes worse when we deal with modulated structures 
with large modulation period. We found that an optimal 
size for the present algorithm is about 200 x 100. The 
total number of Monte Carlo samplings for the given pa- 
rameters of the model a, d and temperature t is about 
10^" 10^^. Parameters of "local heating" are varied in a 
region 2 < Zi,^2 ^ 4. We have checked that a two-times 
increase of the slab size in the x or z-direction does almost 
not change the exponents of power or exponential decay- 
ing or the period of the modulation within the numerical 
accuracy. The result does not depend on the initial config- 
uration of the system. Thus, we believe that the presented 
algorithm allows us to perform accurate numerical simu- 
lations. 



3 Results and their analysis 

We present here detailed results for two values of the 
parameter d dX a = h: d — —0.45 and d — —0.35. We 
choose these values of the parameter d since at least the 
mean-field approximation predicts various ground states 
for zero temperature (commensurate phase with = 6 
for d = —0.45 and ferro-phase for d = —0.35, see Fig. lb). 
The parameter value a = 5 of the DIFFOUR model Q is 
used since the amplitude of modulation then has a value 
comparable with the value of C(0). It makes the analysis 
of the data easier. We have found that the modulation am- 
plitude is much less than the value C(0) for larger values 
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of a (for example, at a > 20). The mean- field approach 
for the DIFFOUR model confirms these observations. 

The correlation functions Cz and are calculated 
at several values of the temperature. The results for Cz 
are presented in Fig. 2 and Fig. 3 for d = —0.45 and 
d = —0.35, respectively. Typical dependencies of Cx for 
both values of the parameter d are shown in Fig. 4. We 
use a log-log scale for the Cx functions. The most infor- 
mative dependencies here are Cz (Fig. 2,3). Qualitatively 
the behaviour of the functions Cx and Cz is different: non- 
monotonic and monotonic dependence on the distance, re- 
spectively. We believe that the influence of the finite-size 
effects starts at a; > 100 and at z > 70 and we do not 
consider these regions in the analysis of the data. The er- 
ror bar for Cx and Cz at moderate values of temperature 
for d = —0.45 and d = —0.35 is about a few percent and 
becomes larger at lower temperatures. 

We analyze the obtained results by fitting the data for 
Cz (Fig. 2,3) using one of the functions H1I2(I . We also try 
to find a power law decay for Cx (Fig. 4) if it is possible. 



3.1 Case d = -0.45 

We start with the case d = —0.45 (Fig. 2). We look for the 
possibility of a phase transition from the para-phase to the 
FIC-phase. The dependencies of the correlation functions 
Cz sA. t > 1.2 can be fitted by equation J^l, i-e. by non- 
monotonic exponential decay. There is no power law decay 
for Cx either for these values of temperature. Thus, we 
conclude that it is floating- fluid phase at t > 1.2 and d — 
-0.45. 

One can fit the dependencies Cz in the vicinity of the 
t ^ 1.2 by nonmonotonic power law decay The results 
of this fitting are presented in Fig. 2 by solid lines. It is pos- 
sible also to find a power law decay for Cx- This region is 
believed to be close to the possible point of the phase tran- 
sition from the floating-fluid phase to the FIC-phase. The 
value of the index ry is approximately equal to 0.25 at the 
temperature ti « 1.2. This coincides with the condition 
of the Kosterlitz-Thouless phase transition. One can find, 
that data for Cz (Fig. 2) for the region 0.55 < < < 1.2 can 
also be fitted by a non-monotonic power law dependence 
|(2Jl with index i] depending on temperature. We plot the 
dependence of the (0.25 — ri{t))'^ on ti — t for the t < 1.2 
(Fig. 5a). One sees a linear dependence, which is a feature 
of the Kosterlitz-Thouless type behaviour ^E]. Thus, the 
phase transition from the floating-fluid phase to the FIC- 
phase at ii « 1.2 is believed to be of Kosterlitz-Thouless 
type. The data for Cx at 0.55 < t < 1.2 can also be fit- 
ted by a monotonic power law decay. Typical dependence 
of Cx for this temperature region is presented in Fig. 4 
{t = 1.0). In our calculations, the indices for Cx and 
Cz do not coincide (the difference is about 10 — 20%). Al- 
though, we don't know whether this difference is of physi- 
cal or numerical origin, we rather attribute it to finite-size 
effects in simulations. Nevertheless, we use here the index 
rj obtained from the data of fitting in the z direction, i.e. 
by equation jSJ for analysis of the data. 



We also show in Fig. 5b how the wavevector changes 
as a function of temperature. We find that the modula- 
tion period is equal to TV = 6 at t < 0.55. The appearance 
of the domain walls (discommensurations - small regions 
with pseudo-period N ^ 6) can be seen at t > 0.55. It 
prevents the system to get a commensurate modulation 
period {N = 6). The point t2 ~ 0.55 is supposed to be 
the point of the phase transition from the FIC-phase to 
the commensurate phase with period = 6 (see Fig. 
5b) . This picture is similar to that for the two-dimensional 
ANNNI model 9 for the transition from the FIC-phase 
to the < 2 > phase. One should note that the accuracy 
of our results is not sufficient to determine whether the 
phase transition from FIC-phase to the < 2 > phase is 
of Pokrovsky-Talapov type or not (Fig. 5). Present re- 
sults are in agreement with previous conclusions by Schulz 
and Haldane et al that only commensurate phase with 
sufficiently small periodicity can transform into a FIC- 
phase 20,21 . It is also possible that narrow commensu- 
rate regions exist in the temperature range 0.55 < t < 1.2. 
Study of this problem demands more accurate methods of 
numerical simulations and analysis of the data. 

3.2 Case d = -0.35 

Before presenting the numerical results, let us stress sev- 
eral features of the d = —0.35 situation. The mean-field 
predicts two transitions: from the para-phase to an incom- 
mensurate one, and than to the ferro-phase. The three- 
dimensional DIFFOUR model indeed behaves so at cer- 
tain parameters On the other hand, the ANNNI model 
does not show two transitions; there is only a single transi- 
tion from the floating-fiuid (disordered) to the ferro-phase. 
For the present system, because of the relatively small 
value a = 5, the modulation is much more pronounced 
and falls off slower than in the ANNNI model. As well, 
the mean-field region for the incommensurate phase is 
quite wide. Therefore one could expect a transition from a 
floating-fiuid to an FIC-phase for d — —0.35. Surprisingly, 
we did not found evidence for this; qualitatively the result 
coincides with that of the ANNNI model. 

Results for the case d = —0.35 for Cz are presented in 
Fig. 3. Data for t > 1.0 can be fltted by non-monotonic 
exponential decay . Results of fltting Cz are shown in 
Fig. 3 by dashed lines. We have performed a detailed stud- 
ies of the region 1.0 < t < 0.8 (part of them is shown in 
Fig. 3). Since the period of the modulation is more than 
14 in this region it is difficult to distinguish whether there 
is a non-monotonic power law decay or not. Besides, the 
fluctuations in this region are increased and our algorithm 
does not allow us to judge the type of decay of the correla- 
tions. We try to fit these dependencies by exponential and 
power law decay. The results for t = 0.95 are presented in 
Fig. 3 by dashed and solid lines, respectively. Both depen- 
dencies show agreement with Monte Carlo results within 
the numerical error. The index i] for power law decay is 
found to be about 0.5. Thus, we believe that there is no 
Kosterlitz-Thouless phase transition for d = —0.35 since 
there is no agreement with a power law decay |2l with the 
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Fig. 2. The correlation function Cz for various temperatures: t — 1.2, t = 1.0, t — 0.8 and t — 0.6. Parameters of the DIFFOUR 
model are a = 5, d — —0.45. Dots - Monte Carlo results, solid lines - results of fitting by non-monotonic power law decay 
The Kosterlitz-Thouless like phase transition takes place at ti ~ 1.2 (index rj = 0.25). 



exponent value rj = 0.25 in this temperature region. One 
can see that the period of the modulation is increased, 
while the amplitude of modulation is decreased. This sit- 
uation is opposite to the case d = — 0.45. A typical result 
for the Cx is presented in Fig. i {t — 0.95). It is difficult 
to find power law decay in the region 1.0 < t < 0.8 for Cx- 
Certainly, there is no index rj = 0.25 in this case either. It 
is difficult to obtain clear dependence of correlations decay 
in the vicinity of the phase transition point due to increase 
of fluctuations and strong influence of the finite-size effects 
(see Fig. 3, t — 0.925). Thus, we conjecture that there is 
only one second order phase transition for d = —0.35 from 
the floating-fluid phase into the ferro-phase. We estimate 
the temperature of this phase transition to be tc ~ 0.91. 
We would like to note that one can suggest a transition 
from floating-fluid phase to a ferro-phase via a phase with 
monotonic exponential decay of correlations with q = 
^ and the existence of a disorder line by analogy with 
the ANNNI model [7j . This phase may exist in a very nar- 
row temperature region at given parameters of the model 
and data for t = 0.925 confirm this. The question of the 
presence of a disorder line in the phase diagram of the two- 
dimensional DIFFOUR model needs an additional study. 

4 Conclusions 

In conclusion, we have studied the DIFFOUR model in 
two dimensions by Monte Carlo simulations. The classi- 
cal DIFFOUR model has two parameters a and d, where 
the latter one is the effective constant of the competition 



of the interactions. Two sets of parameters of the model 
have been considered here in detail: a — 5,d — —0.45 and 
a = 5,d = —0.35. A "local heating" algorithm has been 
used for numerical simulations. This algorithm uses mod- 
ifications of the constants and parameters of the potential 
energy of the DIFFOUR model. Correlation functions in 
both directions have been calculated for analysis of the 
behaviour of the model. 

We have found that correlation functions for the a — 
5,d — —0.45 can be fitted by a non-monotonic expo- 
nential decay at t > 1.2 and by non-monotonic power 
law decay at temperatures 0.55 < t < 1.2. The phase 
transition at ti « 1.2 of Kosterlitz-Thouless type is be- 
lieved to take place with ij — 0.25 from the floating-fluid 
phase to the FIC-phase. The dependence of the index 
ri{t) below the point of the phase transition has a square 
root of ii — t behaviour in agreement with the Kosterlitz- 
Thouless theory. The ground state of the model a.tt < 0.55 
(a = 5, d = —0.45) is a commensurate phase with modu- 
lation period = 6. The second phase transition at these 
parameters of the model from FIC-phase to the commen- 
surate one occurs at t2 ~ 0.55. 

The correlation functions for the case a — 5,d — —0.35 
can be fitted by a non-monotonic exponential decay for 
t > 0.95. Dependencies at i ~ 0.95 can be fitted by non- 
monotonic power law decay also but with 77 > 0.25. Nu- 
merical data do not allow us to distinguish exponential 
from power law decay for this case since the period of mod- 
ulation is too large. Nevertheless, we believe that there 
is no FIC-phase in this case and only one phase transi- 
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Fig. 3. The correlation function Cz for various temperatures: t = 1.2, t — 1.1, t = 1.0, t — 0.95, t = 0.925 and t — 0.8. 
Parameters of the DIFFOUR model are a = 5, d = —0.35. Dots - Monte Carlo results, dashed lines - results of fitting by 
non-monotonic exponential decay (0. Solid line - result of the fitting at t = 0.95 by non-monotonic power law decay index 
rj is approximately equal to 0.5. Single phase transition from floating-fluid phase to the ferro-phase takes place at tc ~ 0.91. 
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Fig. 4. Typical dependencies of the correlation function Cx for d = —0.45 and d = —0.35 (temperatures are equal to t — 1.0 
and t = 0.95, respectively). Logarithmic scale is used for both axes. Results are obtained by Monte Carlo simulations. 
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Fig. 5. a) Dependence of the function (0.25 — on the ti — t below the point of the phase transition from floating-fluid 

phase to the FIC-phase (parameters of the model a = 5, d = —0.45). Dots are the results of the analysis of the Monte Carlo 

data (Fig. 2), dashed line is a linear approximation in the vicinity of the phase transition, b) Dependence of the wavcvoctor q 
on the temperature t (parameters of the model a = 5, d = —0.45). Dots are the results of the analysis of the Monte Carlo data 
(Fig. 2). ti and t-z are the temperatures of the phase transitions from the floating-fluid phase to the FIC-phase and from the 
FIC-phase to the commensurate phase with period N = 6, respectively. 



tion from the floating-fluid phase to the ferro-phase, at 
tc « 0.91. The ground state of the model for f < 0.9 is a 
ferro-phase. 

We would like to point out that the results of the 

Monte Carlo simulations (Fig. 2,3) and the mean-field 
data (Fig. 1) are different. Mean-field approximation pre- 
dicts an incommensurate phase, while Monte Carlo re- 
sults show a FIC-phase with power law decay of corre- 
lations. Mean-fleld theory predicts two phase transitions 
for d — —0.35 and a region of the incommensurate phase, 
while Monte Carlo simulations show only one phase transi- 
tion from the floating-fluid phase to the ferro-phase. Nev- 
ertheless, the low-temperature phase (ferro-phase for d = 
—0.35 and commensurate phase with period = 6 for 
d = —0.45) is predicted correctly in the mean-field ap- 
proximation, as in the limit t = fluctuations are absent 
and the mean- field result is exact. 
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